lavaanパッケージは、Rで共分散構造分析(structual equation modeling)を実行するパッケージです。公式サイトは以下の通りです:
公式サイトにある確認的因子分析とSEMを実行すると、以下のとおりとなる:
library(lavaan)## This is lavaan 0.5-20## lavaan is BETA software! Please report any bugs.xx <- HolzingerSwineford1939
HS_model <- '
visual =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed =~ x7 + x8 + x9
'
fit <- cfa(HS_model, data=xx)
summary(fit)## lavaan (0.5-20) converged normally after 35 iterations
##
## Number of observations 301
##
## Estimator ML
## Minimum Function Test Statistic 85.306
## Degrees of freedom 24
## P-value (Chi-square) 0.000
##
## Parameter Estimates:
##
## Information Expected
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err Z-value P(>|z|)
## visual =~
## x1 1.000
## x2 0.554 0.100 5.554 0.000
## x3 0.729 0.109 6.685 0.000
## textual =~
## x4 1.000
## x5 1.113 0.065 17.014 0.000
## x6 0.926 0.055 16.703 0.000
## speed =~
## x7 1.000
## x8 1.180 0.165 7.152 0.000
## x9 1.082 0.151 7.155 0.000
##
## Covariances:
## Estimate Std.Err Z-value P(>|z|)
## visual ~~
## textual 0.408 0.074 5.552 0.000
## speed 0.262 0.056 4.660 0.000
## textual ~~
## speed 0.173 0.049 3.518 0.000
##
## Variances:
## Estimate Std.Err Z-value P(>|z|)
## x1 0.549 0.114 4.833 0.000
## x2 1.134 0.102 11.146 0.000
## x3 0.844 0.091 9.317 0.000
## x4 0.371 0.048 7.779 0.000
## x5 0.446 0.058 7.642 0.000
## x6 0.356 0.043 8.277 0.000
## x7 0.799 0.081 9.823 0.000
## x8 0.488 0.074 6.573 0.000
## x9 0.566 0.071 8.003 0.000
## visual 0.809 0.145 5.564 0.000
## textual 0.979 0.112 8.737 0.000
## speed 0.384 0.086 4.451 0.000model2 <- '
# measurement model
ind60 =~ x1 + x2 + x3
dem60 =~ y1 + y2 + y3 + y4
dem65 =~ y5 + y6 + y7 + y8
# regressions
dem60 ~ ind60
dem65 ~ ind60 + dem60
# residual correlations
y1 ~~ y5
y2 ~~ y4 + y6
y3 ~~ y7
y4 ~~ y8
y6 ~~ y8
'
fit2 <- sem(model2, data=PoliticalDemocracy)
summary(fit2)## lavaan (0.5-20) converged normally after 68 iterations
##
## Number of observations 75
##
## Estimator ML
## Minimum Function Test Statistic 38.125
## Degrees of freedom 35
## P-value (Chi-square) 0.329
##
## Parameter Estimates:
##
## Information Expected
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err Z-value P(>|z|)
## ind60 =~
## x1 1.000
## x2 2.180 0.139 15.742 0.000
## x3 1.819 0.152 11.967 0.000
## dem60 =~
## y1 1.000
## y2 1.257 0.182 6.889 0.000
## y3 1.058 0.151 6.987 0.000
## y4 1.265 0.145 8.722 0.000
## dem65 =~
## y5 1.000
## y6 1.186 0.169 7.024 0.000
## y7 1.280 0.160 8.002 0.000
## y8 1.266 0.158 8.007 0.000
##
## Regressions:
## Estimate Std.Err Z-value P(>|z|)
## dem60 ~
## ind60 1.483 0.399 3.715 0.000
## dem65 ~
## ind60 0.572 0.221 2.586 0.010
## dem60 0.837 0.098 8.514 0.000
##
## Covariances:
## Estimate Std.Err Z-value P(>|z|)
## y1 ~~
## y5 0.624 0.358 1.741 0.082
## y2 ~~
## y4 1.313 0.702 1.871 0.061
## y6 2.153 0.734 2.934 0.003
## y3 ~~
## y7 0.795 0.608 1.308 0.191
## y4 ~~
## y8 0.348 0.442 0.787 0.431
## y6 ~~
## y8 1.356 0.568 2.386 0.017
##
## Variances:
## Estimate Std.Err Z-value P(>|z|)
## x1 0.082 0.019 4.184 0.000
## x2 0.120 0.070 1.718 0.086
## x3 0.467 0.090 5.177 0.000
## y1 1.891 0.444 4.256 0.000
## y2 7.373 1.374 5.366 0.000
## y3 5.067 0.952 5.324 0.000
## y4 3.148 0.739 4.261 0.000
## y5 2.351 0.480 4.895 0.000
## y6 4.954 0.914 5.419 0.000
## y7 3.431 0.713 4.814 0.000
## y8 3.254 0.695 4.685 0.000
## ind60 0.448 0.087 5.173 0.000
## dem60 3.956 0.921 4.295 0.000
## dem65 0.172 0.215 0.803 0.422SEMでの結果は、多くの場合ダイアグラムで図示することが多いけれど、このlavaanパッケージはplotを自動的には出力しません。
Rでダイアグラムを描くパッケージです。公式サイトは以下のとおりです:
公式サイトにあるサンプルコードを実行すると、以下のとおりです:
library(DiagrammeR)
# Create a node data frame
nodes <-
create_nodes(
nodes = c("a", "b", "c", "d"),
label = FALSE,
type = "lower",
style = "filled",
color = "aqua",
shape = c("circle", "circle",
"rectangle", "rectangle"),
data = c(3.5, 2.6, 9.4, 2.7)
)
edges <-
create_edges(
from = c("a", "b", "c"),
to = c("d", "c", "a"),
rel = "leading_to"
)
graph <-
create_graph(
nodes_df = nodes,
edges_df = edges,
node_attrs = "fontname = Helvetica",
edge_attrs = c("color = blue",
"arrowsize = 2")
)
render_graph(graph)この2つを利用し、lavaanの出力オブジェクトを自動的にDiagrammeRのgraphオブジェクトを生成する関数を作成しました。
create_graph_lavaan <- function(fit, ...) {
library(Hmisc) #順番注意
library(DiagrammeR)
library(dplyr)
# 必要データ取り出し
d <- data.frame(fit@ParTable) %>%
select(lhs,op,rhs,est)
#こっちもありかも
#こっちを使うなら、推定値だけじゃなくてp値も持ってこれるので、
#アスタリスクを自動で付けたりとかできるようになる
#library(lavaan)
#d <- fit %>% parameterEstimates() %>%
# select(lhs,op,rhs,est)
# make NDF
latent <- d %>%
filter(op == "=~") %>%
select(nodes = lhs) %>%
distinct %>%
mutate(shape = "circle", rank = "same")
observed <- d %>%
filter(op != "~1", lhs %nin% latent$nodes) %>%
select(nodes = lhs) %>%
distinct %>%
mutate(shape = "square")
nodes_df <- combine_nodes(latent, observed)
# make EDF
all_paths <- d %>%
filter(op != "~1") %>%
mutate(label = format(round(est, 2),nsmall=2)) %>%
select(-est)
loadings <- all_paths %>%
filter(op == "=~") %>%
mutate(from = lhs, to = rhs, style = "dashed") %>%
select(from, to, style, label)
regressions <- all_paths %>%
filter(op == "~~") %>%
mutate(to = lhs, from = rhs, style = "solid", dir = "both") %>%
select(from, to, style, label, dir)
predictions <- all_paths %>%
filter(op == "~") %>%
mutate(from = lhs, to = rhs, style = "solid") %>%
select(from, to, style)
edges_df <- combine_edges(loadings, regressions, predictions)
graph <- create_graph(
nodes_df = nodes_df,
edges_df = edges_df,
graph_attrs = c("layout = 'dot'")
)
return(graph)
}先のlavaanで実行したfitとfit2を実際に実行すると、以下のとおりです:
DgmR1 <- create_graph_lavaan(fit)## Loading required package: lattice## Loading required package: survival## Loading required package: Formula## Loading required package: ggplot2##
## Attaching package: 'Hmisc'## 以下のオブジェクトは 'package:base' からマスクされています:
##
## format.pval, round.POSIXt, trunc.POSIXt, units##
## Attaching package: 'dplyr'## 以下のオブジェクトは 'package:Hmisc' からマスクされています:
##
## combine, src, summarize## 以下のオブジェクトは 'package:stats' からマスクされています:
##
## filter, lag## 以下のオブジェクトは 'package:base' からマスクされています:
##
## intersect, setdiff, setequal, unionDgmR2 <- create_graph_lavaan(fit2)
render_graph(DgmR1)render_graph(DgmR2)現在の状態では配置に問題があるため、改善が必要となります。またCFAと普通のSEMにしか対応していないので、他母集団SEMに対しても対応させたい。
以上です。